clear all;clc
addpath(genpath('.'));  %添加工具箱路径
load R_Scalar
load Scale
load id
load step
%% 超参数设置
% id=2;  
% bitstream=[];%若干行，两列，第一列代表水印的提取信息，第二列代表这个水印的置信度
% Bit_Stream=[];
if id==9
    id=id-1;
end
if id==1
    id=id+1;
end
%% 获取最佳尺度的前后两个参考Map
A=R_Scalar{id};%最佳置信Map
B=R_Scalar{id+1}; %辅助参考Map B
C=R_Scalar{id-1}; %辅助参考Map C
% A=B;
SID=[id,id+1,id-1];
Scalar={A,B,C};
count=1;Intensity=[];%计算内部区域面积
for data=Scalar
    A=data{1};
%% 可能需要做个阈值化，提前去除一些异常值,还是得在这一步做一个中心内部去异常点，或者在提取时做掩膜

%% 取前14强波峰和前14强波谷，这个可以暴力遍历一下
    p_count=1;
for max_peak=13:18 %这个依旧是个头痛的问题13.18
    peak_loc=[];
    valley_loc=[];
    A=data{1};%重载A
    for i=1:max_peak%大不了这个也遍历一下呗，找个面积最大的即可，目前测试的最优解范围13~18
        %波峰
        [~,mpos]=max(A(:));[row,col] = ind2sub(size(A),mpos);
        peak_loc=[peak_loc;[row,col]];
        [row,col]=check_boundary(row,col,A);
        A(row-3:row+3,col-3:col+3)=0;%这个地方得改一下，避免越界
        %波谷
        [~,mpos]=min(A(:));[row,col] = ind2sub(size(A),mpos);
        valley_loc=[valley_loc;[row,col]];
        [row,col]=check_boundary(row,col,A);
        A(row-3:row+3,col-3:col+3)=0;
    end

%% 均值分割四条线,需要进行离散值清理，but还没做，要么在前做，要么在后做
sorted_peak=sortrows(peak_loc,2);sorted_valley=sortrows(valley_loc,1); 
orgin_sorted_peak=sorted_peak;orgin_sorted_valley=sorted_valley;
Line_A1=[];Line_A2=[];Line_B1=[];Line_B2=[];
P_m=mean(sorted_peak(:,2));V_m=mean(sorted_valley(:,1));
sorted_peak(sorted_peak(:,2)<P_m)=0;temp=find(sorted_peak~=0);ini_row=temp(1);
sorted_valley(sorted_valley(:,1)<V_m)=0;temp=find(sorted_valley~=0);ini_col=temp(1);
Line_A1=orgin_sorted_peak(1:ini_row-1,:);
Line_A2=orgin_sorted_peak(ini_row:end,:);
Line_B1=orgin_sorted_valley(1:ini_col-1,:);
Line_B2=orgin_sorted_valley(ini_col:end,:);

%% 依次线性拟合四条线
% a1=[k,b]
%% 采用简单线性回归，只要斜率和截距即可
a1=polyfit(Line_A1(:,1),Line_A1(:,2),1);%1阶多项式最小二乘拟合

% y1=polyval(a1,Line_A1(:,1));
% plot(Line_A1(:,1),Line_A1(:,2),'*',Line_A1(:,1),y1,'g-','linewidth',2);
a2=polyfit(Line_A2(:,1),Line_A2(:,2),1);%1阶多项式最小二乘拟合
b1=polyfit(Line_B1(:,1),Line_B1(:,2),1);%1阶多项式最小二乘拟合
b2=polyfit(Line_B2(:,1),Line_B2(:,2),1);%1阶多项式最小二乘拟合
%% 采用带误差估计的线性回归
% [a1,s]=polyfit(Line_A1(:,1),Line_A1(:,2),1);%1阶多项式最小二乘拟合 具有误差估计
% [a1,delta] = polyval(a1,Line_A1(:,1),s);
% [a2,s]=polyfit(Line_A2(:,1),Line_A2(:,2),1);%1阶多项式最小二乘拟合 具有误差估计
% [a2,delta] = polyval(a2,Line_A2(:,1),s);
% [b1,s]=polyfit(Line_B1(:,1),Line_B1(:,2),1);%1阶多项式最小二乘拟合 具有误差估计
% [b1,delta] = polyval(b1,Line_B1(:,1),s);
% [b2,s]=polyfit(Line_B2(:,1),Line_B2(:,2),1,1);%1阶多项式最小二乘拟合 具有误差估计
% [b2,delta] = polyval(b2,Line_B2(:,1),s);
%% 求四个顶点
[a1_b1_x,a1_b1_y]=linecross(a1(1),a1(2),b1(1),b1(2));%左上
[a1_b2_x,a1_b2_y]=linecross(a1(1),a1(2),b2(1),b2(2));%左下
[a2_b1_x,a2_b1_y]=linecross(a2(1),a2(2),b1(1),b1(2));%右上
[a2_b2_x,a2_b2_y]=linecross(a2(1),a2(2),b2(1),b2(2));%右下
Line_Set{count,p_count}={[a1_b1_x,a1_b1_y],[a1_b2_x,a1_b2_y],[a2_b1_x,a2_b1_y],[a2_b2_x,a2_b2_y]};
if check_Line([[a1_b1_x,a1_b1_y];[a1_b2_x,a1_b2_y];[a2_b1_x,a2_b1_y];[a2_b2_x,a2_b2_y]],A)
    %Intensity利用多边形面积&夹角比例(作为惩罚项)
    sub_pix_location=[[a1_b1_x,a1_b1_y];[a1_b2_x,a1_b2_y];[a2_b1_x,a2_b1_y];[a2_b2_x,a2_b2_y]];
    subrow=sub_pix_location(:,1);
    subcol=sub_pix_location(:,2);
    row_2=[subrow(1),subrow(3),subrow(4),subrow(2),subrow(1)];%顺时针
    col_2=[subcol(1),subcol(3),subcol(4),subcol(2),subcol(1)];
    Angle_Punishment=cac_angle_intensity([row_2,subrow(3)],[col_2,subcol(3)]);
%     Angle_Punishment=1;
    Intensity(count,p_count)=polyarea(col_2,row_2)/(Scale(SID(count)))*Angle_Punishment;
%     Intensity(count,p_count)=norm([a1_b1_x,a1_b1_y]-[a2_b2_x,a2_b2_y])*norm([a1_b2_x,a1_b2_y]-[a2_b1_x,a2_b1_y])/(Scale(SID(count)));
else
    Intensity(count,p_count)=0;
end
p_count=p_count+1;
end
count=count+1;
end
[~,mpos]=max(Intensity(:));[p_xid,p_yid] = ind2sub(size(Intensity),mpos);
Point_Set=Line_Set{p_xid,p_yid};
save Point_Set Point_Set
save SID SID
save p_xid p_xid
save p_yid p_yid
%看看效果
figure,imshow(R_Scalar{SID(p_xid)},[]);
hold on;
% plot(a1_b1_y,a1_b1_x,'o','LineWidth',8);
% plot(a1_b2_y,a1_b2_x,'o','LineWidth',8);
% plot(a2_b1_y,a2_b1_x,'o','LineWidth',8);
% plot(a2_b2_y,a2_b2_x,'o','LineWidth',8);
pix_location=[];
for location=Point_Set
    location=location{1};
    pix_location=[pix_location;(location+1)*step];
    plot(location(2),location(1),'o','LineWidth',8);
end
save pix_location pix_location
% Stego=rgb2ycbcr(imread('./Phone_image/loc_3.jpg'));
% Stego=Stego(:,:,1);
% left_up=[(a1_b1_x-1)*64+1,(a1_b1_y-1)*64+1];
% right_up=[(a1_b2_x-1)*64+1,(a1_b2_y-1)*64+1];
% left_down=[(a2_b1_x-1)*64+1,(a2_b1_y-1)*64+1];
% right_down=[(a2_b2_x-1)*64+1,(a2_b2_y-1)*64+1];
% figure,imshow(Stego,[]);
% hold on;  %光plot是不行的，必须要有这句，才能实现在 原图 标记点这一目标
% plot(left_up(2),left_up(1),'o');  %用来标记点的核心语句
% plot(right_up(2),right_up(1),'o');%红蓝圈就是标记的点。颜色没设置会随机分配
% plot(left_down(2),left_down(1),'o');  %用来标记点的核心语句
% plot(right_down(2),right_down(1),'o');%红蓝圈就是标记的点。颜色没设置会随机分配
% point=[left_up;right_up;left_down;right_down];
% 看看效果
% line([a1_b1_x,a2_b1_x],[a1_b1_y,a2_b1_y]);%line([x1,x2],[y1,y2])
% hold on
% line([a1_b1_x,a2_b2_x],[a1_b1_y,a2_b2_y]);
% hold on
% line([a1_b2_x,a2_b1_x],[a1_b2_y,a2_b1_y]);
% hold on
% line([a1_b2_x,a2_b2_x],[a1_b2_y,a2_b2_y]);
% hold on
% line([a1_b1_x,a1_b2_x],[a1_b1_y,a1_b2_y]);
% hold on
% line([a2_b1_x,a2_b2_x],[a2_b1_y,a2_b2_y]);
% hold off
% Line_A=
% for i=1:15
%     if sorted_peak(i,1)<P_m
%         
%     else
%         
%     end
%     if sorted_valley(i,1)<V_m
%         
%     else
%         
%     end
% end
% %% 均值方差分割图像
% Mean=mean2(A);
% Var=std2(A);
% location_1=A>(Mean+Var);
% location_2=A<(Mean-Var);
% Location_ZOOM=location_1+location_2;
% Location_ZOOM(Location_ZOOM>1)=1;
% %% 可能需要做一个霍夫变换校正水印图像，当拍的很歪的时候
% trans_detetion_img = bwmorph(bwmorph(Location_ZOOM,'dilate',2),'erode',2);  %膨胀
% % figure,
% % subplot(121),imshow(Location_ZOOM,[]); title('无形态学')
% % subplot(122),imshow(trans_detetion_img,[]);title('形态学')
% %% 获取目标水印块区域
% row_num=sum(Location_ZOOM,2);row_num(row_num<15)=0;%每行和,这里的10代表一排最多是个波峰波谷算是个先验
% col_num=sum(Location_ZOOM,1);col_num(col_num<15)=0;%每列和
% temp=find(row_num~=0);ini_row=temp(1);end_row=temp(end);
% temp=find(col_num~=0);ini_col=temp(1);end_col=temp(end);
% detect_region=imcrop(A,[ini_col,ini_row,end_col-ini_col-0,end_row-ini_row-0]);
% orgin_detect_region=detect_region;
% %% 获取原始水印长度，依据row_num & col_num
% [~,LOCS_row]=findpeaks(row_num);[~,LOCS_col]=findpeaks(col_num);
% bit_row=get_bitlength(LOCS_row);bit_col=get_bitlength(LOCS_col);

%% 横刀切面法，问题显著
% Bit_map=zeros(size(detect_region));
% while (1)
% [Max_layer]=max(max(detect_region));
% [Min_layer]=min(min(detect_region));
% % 领域搜索，更新detect_region，将以max_id,min_id 的3*3领域内的值记为0
%     if Max_layer>0
%         Bit_map(Max_layer==detect_region)=1;
%         [row,col]=find(Max_layer==detect_region);
%         detect_region(row-3:row+3,col-3:col+3)=0;
%     end
%     if Min_layer<0
%         Bit_map(Min_layer==detect_region)=-1;
%         [row,col]=find(Min_layer==detect_region);
%         detect_region(row-3:row+3,col-3:col+3)=0;
%     end
% end
%% 尝试利用一个4*4的窗口来提取水印信息，也可能3*3，自己设计的东西过多，泛化性能可能不太好
% [row,col]=size(detect_region);
% detect_region=My_padding(detect_region,windows_size);
%% 二值化,有必要吗？有必要，但是只需要一值化，就是将梯度小的地方置零，其次再将block少的地方置零
% Var_region=std2(detect_region)/3;Mean_region=mean2(detect_region);
% % Location_1=detect_region>(Mean_region+Var_region);%代表1
% % Location_2=detect_region<(Mean_region-Var_region);%代表0
% Location_0=detect_region>(Mean_region-Var_region) & detect_region<(Mean_region+Var_region);%代表null
% % detect_region(Location_1)=1;detect_region(Location_2)=-1;
% detect_region(Location_0)=0;
% C_row=sum(detect_region==0,2);
% C_col=sum(detect_region==0,1);
% % [5,9]
% [~,LOCS_row]=findpeaks(C_row,'minpeakheight',mean(C_row));[~,LOCS_col]=findpeaks(C_col,'minpeakheight',mean(C_col));
% %% 极值检查，也可以用来筛选尺度
% detect_region(LOCS_row,:)=0;detect_region(:,LOCS_col)=0;
% 
% %% 水印提取
% Bit_Stream=extract_msg(orgin_detect_region,LOCS_row,LOCS_col);
%  Bit_Stream=Bit_Stream';
% end
% for i=1:row
%     for j=1:col
%         judge_mont=detect_region((i:i-1+windows_size),(j:j-1+windows_size));
%         sum_ele=sum(judge_mont(:));confident=abs(sum_ele/(windows_size^2));
%         if sum_ele >0
%            bitstream=[bitstream,[1,confident]];
%         end
%         if sum_ele <0
%            bitstream=[bitstream,[0,confident]];
%         end
%         if sum_ele ==0
%             bitstream=[bitstream,[-1,0]];
%         end
%     end
%     Bit_Stream=[Bit_Stream;bitstream];
%     bitstream=[];
% end
% Bit_Set=Bit_Stream(:,1:2:end);
% Confident_Set=Bit_Stream(:,2:2:end);
% 
% step_row=size(Bit_Set,1)/bit_row;step_col=size(Bit_Set,2)/bit_col;
% S1=floor(step_row);S2=floor(step_col);
% Rate_S1=S1-step_row;Rate_S2=S2-step_col;
% for i=1:bit_row
%     for j=1:bit_col
%          block_bit=Bit_Set((i-1)*S1+1:i*S1,(j-1)*S2+1:j*S2);
%          block_confident=Confident_Set((i-1)*S1+1:i*S1,(j-1)*S2+1:j*S2);
%          % 统计1 和 -1的数量
%          a=sum(block_bit==1); b=sum(block_bit==0);
%          if a-b>6
%              bit=1;
%          else
%             if b-a>6
%              bit=0;
%             else %考察confident
%                 if block_bit(confident==max(max(confident)))>0
%                     bit=1;
%                 else
%                     bit=0;
%                 end
%             end
%          end
%          bitstream=[bitstream,bit];
%     end
% end
%% 尝试利用一维波峰波谷来提取水印，按行波峰法和列波峰法不太行
% [w,h]=size(detect_region);
% for i=1:w
%     detect_row=detect_region(i,:);
%     findpeaks(detect_row);
% end
% Mean_S=mean2(detect_region);Var_S=std2(detect_region);
% location_peak=detect_region>(Mean+Var/4);location_valley=detect_region<(Mean-Var/4);
% % subplot(121),imshow(location_row);subplot(122),imshow(location_col);% 检查一下
% row_num=sum(location_peak,2);row_num(row_num<10)=0;%每行和,这里的10代表一排最多是个波峰波谷
% col_num=sum(location_peak,1);col_num(col_num<10)=0;%每列和
% [~,LOCS_row]=findpeaks(row_num);[~,LOCS_col]=findpeaks(col_num);
% bit_row=get_bitlength(LOCS_row);bit_col=get_bitlength(LOCS_col);


% 生成原始比特，计算Ber
% load key
% load p
% rng(key); %%key为伪随机序列种子
% msg = num2str(uint8(rand(1,100*100)));
% msg=strrep(msg,' ','');%剔除空格
% embed_msg=msg(1:p-1);%获取嵌入的秘密信息
% embed_msg=mychar2double(embed_msg);
% Final_bit=1-Bit_Stream(:);
% diff_bit_2=sum(abs(Final_bit-embed_msg));
% NC=1-diff_bit_2/length(embed_msg)
% % 计算Ber
% bitstream=bitstream';
% embed_msg=mychar2double(embed_msg);
% diff_bit_2=sum(abs(bitstream-embed_msg));
% NC=1-diff_bit_2/length(embed_msg)
% window_size=4;
% for 
% size()
%% 均值方差分割图像
% Mean=mean2(A);
% Var=std2(A);
% location_1=A>(Mean+Var);
% location_2=A<(Mean-Var);
% Location_ZOOM=location_1+location_2;
% Location_ZOOM(Location_ZOOM>1)=1;
% cv.findContours(Location_ZOOM);
% Trans1=hough_correction(Location_ZOOM);
% Trans2=hough_correction(trans_Location_ZOOM);
% figure,
% subplot(211),imshow(trans_Location_ZOOM);
% subplot(212),imshow(Location_ZOOM);
% figure,
% subplot(211),imshow(A>(Mean+Var));
% subplot(212),imshow(A<(Mean-Var));
% trans_detetion_img = bwmorph(bwmorph(Location_ZOOM,'dilate',2),'erode',2);  %膨胀